function diff = model_SA_inner(theta,param,aux)
    
    param.c=abs(theta);
    param = model_2(param,aux);
    
	m   = linspace( param.m_l, param.m_h, aux.n)';
    q_p = q_p_fun( m, param);
    w = q_p./length(q_p).*(param.m_h-param.m_l);
    J = J_SA( m, param);
    % Normalize vacancy filling rate to one in initial steady state
    c = param.c_v.*q_fun(param.m_u,param)./q_fun(param.m_h,param)./chi_fun(param) ...
        + w'*J./q_fun(param.m_h,param);
    diff = c-param.c;

end